#ifndef DYN_ODE_DOPRI54_H
#define DYN_ODE_DOPRI54_H

#include "dyn_ode.h"

/*! \brief Implementation of the Dormand-Prince
 * method of order 5(4).
 *
 * This method uses a fifth-order
 * Runge-Kutta algorithm along with an
 * embedded fourth-order method. The difference
 * between the two methods gives an error
 * estimate. Based on the error, we can
 * modify the stepsize in order to stay
 * within the desired accuracy.
 */

class DOPRI54 : public ODE
{
public:
    /*! The constructor setups the variables
     * of the integrator.
     */
    DOPRI54(MultiFunctor &func,
            colvec _initCond,
            double _start,
            double _end,
            double _initStepsize,
            double _minStepsize,
            double _maxStepsize,
            double _absTol,
            double _relTol,
            int _maxIterations);

    /*! @name Inherited virtual functions.
     * Implementation of the inherited pure
     * virtual functions.
     */
    //@{
    int integrate();
    colvec next(colvec previousStep);
    //@}

    /*! @name Accessor Functions
     * We access the old time step and the
     * PI control scheme variables.
     */
    //@{
    double getOldStep(){return oldStep;}
    double getBeta(){return beta;}
    double getAlpha(){return alpha;}
    double getMinScale(){return minScale;}
    double getMaxScale(){return maxScale;}

    void setOldStep(double _oldStep){oldStep=_oldStep;}
    void setBeta(double _beta){beta=_beta;}
    void setAlpha(double _alpha){alpha=_alpha;}
    void setMinScale(double _minScale){minScale=_minScale;}
    void setMaxScale(double _maxScale){maxScale=_maxScale;}
    //@}

protected:
    /*! Size of the step we have taken in
     * a single iteration.
     */
    double oldStep;

    /*! @name PI Stepsize Control
     * We hold the last error we have in memory
     * to make the stepsize correction.
     * See Press et al, 2007, Numerical Recipes
     * for details about the initialization
     * values.
     */
    //@{
    double errorOld=1e-4;
    double beta=0.08;
    double alpha=0.2-0.75*beta;
    double minScale=0.2;
    double maxScale=10.0;
    //@}


    /*! @name Butcher Tableau
     * The following constants form the
     * Butcher tableau of the method.
     */
    //@{
    /*! @name Runge-Kutta matrix.*/
    //@{
    const double a21 = 1./5.;
    const double a31 = 3./40.;
    const double a32 = 9./40.;
    const double a41 = 44./45.;
    const double a42 = -56./15.;
    const double a43 = 32./9.;
    const double a51 = 19372./6561.;
    const double a52 = -25360./2187.;
    const double a53 = 64448./6561.;
    const double a54 = -212./729.;
    const double a61 = 9017./3168.;
    const double a62 = -355./33;
    const double a63 = 46732./5247.;
    const double a64 = 49./176.;
    const double a65 = -5103./18656.;
    //@}
    /*! @name Weights (actual value) */
    //@{
    const double b1 = 35./384.;
    const double b2 = 0.0;
    const double b3 = 500./1113.;
    const double b4 = 125./192.;
    const double b5 = -2187./6784.;
    const double b6 = 11./84.;
    //@}
    /*! @name Weights (error estimate) */
    //@{
    const double c1 = 5179./57600.;
    const double c2 = 0.0;
    const double c3 = 7571./16695.;
    const double c4 = 393./640.;
    const double c5 = -92097./339200;
    const double c6 = 187./2100.;
    //@}
    //@}
};

#endif // DYN_ODE_DOPRI54_H
